Single-cell RNA-Seq adult dataset cleanup

Loading the required packages for the rest of the script. Run once per session.

library(Seurat)
library(sctransform)
library(gridExtra)
library(dplyr)
library(tximport)
library(Matrix)
library(ggplot2)
library(stringr)
library(SeuratWrappers)

0.1 Read big Seurat

Here we read the unfiltered Seurat object with 80,235 cells, processed with SCTransform. And the filtered Seurat object with 52,638 cells, processed with SCTransform. The filtered version we used a cutoff of nFeature_RNA > 800 and percent.mt < 15.

For both, UMAp, Neighbors and Clusters were calculated.

0.2 Sample Big Seurat into Small Seurat (10,000 cells, named pleuro.2)

Running FindAllMarkers function from Seurat in order to look at clusters of bad quality on a small subset of cells that are evenly and randomly sampled from the big dataset to avoid memory issues and to optimize the code locally before running it on HPC. You won’t have to do this on the DEVELOPMENTAL dataset. For screening bad clusters, do FeaturePlot(object, features = c(“percent.mt”,“nFeature_RNA”)) and take those that have high percent.mt and low nFeature_RNA to train your classifier.

# You can run this part as an example of how it works
pleuro.2 <- pleuro_filter[, sample(colnames(pleuro_filter), size = 10000, replace=F)]

#FeaturePlot(new_object, cells = names(new_object@active.ident[new_object@active.ident %in% c("17","30")]), label = T,label.size = 1, pt.size = 0.3, reduction = 'umap', features = c("SATB1","SST")) + NoLegend()

FeaturePlot(pleuro_unfilter, label = T, label.size = 4, pt.size = 0.3, reduction = 'umap', features = c("percent.mt","nFeature_RNA"), ncol = 2) + NoLegend()
FeaturePlot(pleuro_filter, label = T, label.size = 4, pt.size = 0.3, reduction = 'umap', features = c("percent.mt","nFeature_RNA"), ncol = 2) + NoLegend()

You can also run FindAllMarkers function to see whether you can identify clusters that are driven by the expression of mitochondrial genes and/or ribosomal genes and include them as “bad” clusters for the classifier. This function usually takes like 8 hrs locally so I would run it on Ginsburg if you are really invested.

all.markers <- FindAllMarkers(pleuro.2)

FeaturePlot(pleuro.2, label = T, label.size = 4, pt.size = 0.3, reduction = 'umap', features = c("percent.mt","nFeature_RNA"), ncol = 2) + NoLegend()

Screening for good and bad cluster on the basis of marker genes. Some clusters do not have good marker genes to be defined as clusters and/or their differences are being driven by ribosomal/mitocondrial genes.

cluster0 <- subset(all.markers, cluster == 0, select = c(gene, pct.1, pct.2)) # Good
cluster0[,4] <- cluster0$pct.1 - cluster0$pct.2

cluster1 <- subset(all.markers, cluster == 1, select = c(gene, pct.1, pct.2)) # Good cluster
cluster1[,4] <- cluster1$pct.1 - cluster1$pct.2

cluster6 <- subset(all.markers, cluster == 6, select = c(gene, pct.1, pct.2)) # 
cluster6[,4] <- cluster6$pct.1 - cluster6$pct.2

cluster17 <- subset(all.markers, cluster == 17, select = c(gene, pct.1, pct.2)) # BAD RIBO
cluster17[,4] <- cluster17$pct.1 - cluster17$pct.2

cluster21 <- subset(all.markers, cluster == 21, select = c(gene, pct.1, pct.2)) # Good
cluster21[,4] <- cluster21$pct.1 - cluster21$pct.2

cluster24 <- subset(all.markers, cluster == 24, select = c(gene, pct.1, pct.2))
cluster24[,4] <- cluster24$pct.1 - cluster24$pct.2

cluster26 <- subset(all.markers, cluster == 26, select = c(gene, pct.1, pct.2))
cluster26[,4] <- cluster26$pct.1 - cluster26$pct.2

cluster29 <- subset(all.markers, cluster == 29, select = c(gene, pct.1, pct.2))
cluster29[,4] <- cluster29$pct.1 - cluster29$pct.2

cluster40 <- subset(all.markers, cluster == 40, select = c(gene, pct.1, pct.2))
cluster40[,4] <- cluster40$pct.1 - cluster40$pct.2

cluster41 <- subset(all.markers, cluster == 41, select = c(gene, pct.1, pct.2))
cluster41[,4] <- cluster41$pct.1 - cluster41$pct.2

cluster50 <- subset(all.markers, cluster == 50, select = c(gene, pct.1, pct.2))
cluster50[,4] <- cluster50$pct.1 - cluster50$pct.2

cluster52 <- subset(all.markers, cluster == 52, select = c(gene, pct.1, pct.2))
cluster52[,4] <- cluster52$pct.1 - cluster52$pct.2

0.3 Creating meta.data for the classifier

I will use all the bad clusters from this small subset of cells (10,000 cells). For the method to work properly, you have to select all the clusters that look bad, since it will consider as “good” all the rest of the cells. Once you have identified all the bad clusters, make sure to write all their identities in the code. You will want to train your set with ~10% of the total number of “bad cells” that you have (look at length(bad.names)) or table(object@meta.data$seurat_clusters) and sum the bad clusters).

train.final = data.frame(train, ytrain=as.factor(ytrain))

Error in data.frame(train, ytrain = as.factor(ytrain)) : 
  arguments imply differing number of rows: 700, 11000

0.4 Training the classifier

library(e1071)
svmfit=svm(ytrain~., data=train.final, kernel="linear", cost=0.1, scale=F)

set.seed (1)
tune.out=tune(svm,ytrain~.,data=train.final,kernel="linear", ranges=list(cost=c(0.1,1,10,100,500)))
summary(tune.out)


bestmod=tune.out$best.model
summary(bestmod)

test.bad.names<-bad.names[!bad.names %in% train.bad.names] 
test.good.names<-good.names[!good.names %in% train.good.names] 
test.bad<-data.info[test.bad.names,] 
test.good<-data.info[test.good.names,] 
test<-rbind(test.good,test.bad) 

test<-data.info[!rownames(data.info) %in% rownames(train),] 

test<-as.matrix(test)
ytest<-matrix(c(rep(1,length(test.good.names)),rep(-1,length(test.bad.names)))) 

test.final = data.frame(test, ytest=as.factor(ytest))

ypred=predict(bestmod,test.final)
table(predict=ypred, truth=test.final$ytest)

ypred<-as.data.frame(ypred)

0.5 Adding the SVM classification to Seurat Object.

pleuro.2$ypred <- ypred

pred<-as.data.frame(pleuro.2$ypred)
rownames(pred)<-rownames(pleuro.2[[]])
pred[train.bad.names,]<-"-1"
pred[train.good.names,]<-"1"
colnames(pred)<-c("svm_class")
pleuro.2$svm_class <- pred

DimPlot(pleuro.2, split.by="svm_class", pt.size=0.7, label = T, group.by = "library") #+ NoLegend()

0.6 How many cells were thrown away?

# saveRDS(pleuro.2, file = "pleuro2_cleanup.rds")
table(pleuro.2@meta.data$svm_class)

  -1    1 
3110 6890 

pleuro.2 <- RunUMAP(pleuro.2, dims=1:ndims)

13:45:53 UMAP embedding parameters a = 0.9922 b = 1.112
13:45:53 Read 6890 rows and found 120 numeric columns
13:45:53 Using Annoy for neighbor search, n_neighbors = 30
13:45:53 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:45:54 Writing NN index file to temp file /var/folders/h4/17k58rr56z56y6dsspm9qmn40000gp/T//RtmpbeyK9U/file76848e7eb26
13:45:54 Searching Annoy index using 1 thread, search_k = 3000
13:45:56 Annoy recall = 100%
13:45:57 Commencing smooth kNN distance calibration using 1 thread
13:45:58 Initializing from normalized Laplacian + noise
13:45:59 Commencing optimization for 500 epochs, with 270056 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:46:09 Optimization finished

pleuro.2 <- RunUMAP(pleuro.2, dims=1:ndims)

13:46:09 UMAP embedding parameters a = 0.9922 b = 1.112
13:46:09 Read 6890 rows and found 120 numeric columns
13:46:09 Using Annoy for neighbor search, n_neighbors = 30
13:46:09 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:46:10 Writing NN index file to temp file /var/folders/h4/17k58rr56z56y6dsspm9qmn40000gp/T//RtmpbeyK9U/file76871ff312e
13:46:10 Searching Annoy index using 1 thread, search_k = 3000
13:46:12 Annoy recall = 100%
13:46:12 Commencing smooth kNN distance calibration using 1 thread
13:46:14 Initializing from normalized Laplacian + noise
13:46:14 Commencing optimization for 500 epochs, with 270056 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
13:46:24 Optimization finished

pleuro.2 <- FindNeighbors(pleuro.2, dims=1:ndims)

Computing nearest neighbor graph

Computing SNN

pleuro.2 <- FindClusters(pleuro.2, resolution = 2)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 6890
Number of edges: 234564

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8634
Number of communities: 43
Elapsed time: 0 seconds

DimPlot(pleuro.2, reduction="umap", label=T, pt.size = 0.3) + NoLegend()

I ran AllMarkers as a form of sanity check, just to see whether cluster 0 (which seemed like a blob in the middle) was a good cluster or not. Do not run this locally.

all.markers <-FindAllMarkers(pleuro.2)
cluster0 <- subset(all.markers, cluster == 0, select = c(gene, pct.1, pct.2)) # Bad cluster
cluster0[,4] <- cluster0$pct.1 - cluster0$pct.2

Now it is a good cluster.

0.7 From here, I used the above code but with the large datasets (80,235 and 52,638)

I used the exact same code, only that this time I used the large Seurat objects and ran everything on Ginsburg. I am reading the output file from that code with 36,116 cells for pre-filtered and 41,657 for unfiltered. For this part, I tried different nPCs (ndims) in the code with different resolution.

After looking at the code together and long deliberation, we chose to stick to cluster resolution = 2. Next, we will subset this object to keep only the neurons. The rest of the work to be done with this object is showing neuronal and non-neuronal cells in a general way.

png(paste0("UMAP_byanimal.png"), width=48, height=12, units = "in", res=400)
DimPlot(pleuro, reduction="umap", label=F, pt.size = 0.1, split.by = "animal") + NoLegend()
dev.off()

null device 
          1 

pdf(paste0("UMAP_byanimal.pdf"), width=48, height=12)
DimPlot(pleuro, reduction="umap", label=F, pt.size = 0.1, split.by = "animal") + NoLegend()
dev.off()

null device 
          1 

1 FindAllMarkers for Subset pleuro (37K)

# RUN IN CLUSTER
all.markers <- FindAllMarkers(pleuro)

2 FeaturePLots

A couple more FeaturePlots as we continue to check whether everything is ok.

3 Quality control of the final dataset

3.1 VlnPlot

Violin plots for nFeature, nCount and percent.mt grouped by animal and by library, in the cleaned dataset.

png(paste0("Vln_quality_cluster2.png"), width=10, height=6, units = "in", res=400)
VlnPlot(pleuro, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), pt.size = 0, group.by = "animal", ncol = 3)
dev.off()

null device 
          1 

pdf(paste0("Vln_quality_cluster2.pdf"), width=10, height=6)
VlnPlot(pleuro, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), pt.size = 0, group.by = "animal", ncol = 3)
dev.off()

null device 
          1 

VlnPlot(pleuro, features = c("nFeature_RNA", "nCount_RNA"), group.by = "library", ncol = 2)

saveRDS(pleuro, file = "pleuro_f_final.rds")
#saveRDS(pleuro.uf, file = "pleuro_uf_final.rds")

Getting the information by cluster.

all.markers <- read.csv(file = "all_markers_neurons_res2.csv")
list_markers <- list()
a <- 0
for(i in 1:max(all.markers$cluster)) {
  cluster <- subset(all.markers, cluster == a, select = c(gene, pct.1, pct.2))
  cluster[,4] <- cluster$pct.1 - cluster$pct.2
  list_markers[[i]] <- cluster
  a <- a + 1
}

4 Generating a merged salamander dataset by general cell type

new.cluster.ids <- c('Excitatory telencephalic' #0
,'Excitatory telencephalic'
,'Excitatory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Excitatory telencephalic'
,'Inhibitory telencephalic' #10
,'Immature neurons'
,'Excitatory telencephalic' #12
,'Inhibitory telencephalic'
,'Excitatory telencephalic'
,'Excitatory telencephalic' 
,'Excitatory telencephalic'
,'Excitatory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic' #19
,'Oligodendrocyte precursor cells' 
,'Ependymoglia'
,'Inhibitory telencephalic' #22
,'Ependymoglia'
,'Inhibitory telencephalic' #24
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Non-telencephalic'
,'Inhibitory telencephalic'
,'Excitatory telencephalic'
,'Excitatory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic' #33
,'Microglia' 
,'Excitatory telencephalic'#35
,'Ependymoglia'
,'Excitatory telencephalic' #37
,'Excitatory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Excitatory telencephalic'
,'Inhibitory telencephalic' #43
,'Vascular leptomeningeal cells' 
,'Inhibitory telencephalic' #45
,'Non-telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Non-telencephalic' 
,'Inhibitory telencephalic'
,'Inhibitory telencephalic' #52
,'Immature neurons'
,'Inhibitory telencephalic' #54
,'Inhibitory telencephalic'
,'Excitatory telencephalic'
,'Excitatory telencephalic'
,'Inhibitory telencephalic'
,'Non-telencephalic'
,'Inhibitory telencephalic' #60
,'Olfactory ensheating cells'
,'Inhibitory telencephalic' #62
,'Perivascular macrophages'
,'Inhibitory telencephalic' #64
,'Inhibitory telencephalic'
,'Oligodendrocytes'
,'Non-telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic'
,'Inhibitory telencephalic')
names(new.cluster.ids) <- levels(pleuro)
pleuro.dotplot <- RenameIdents(pleuro, new.cluster.ids)

#saveRDS(pleuro.dotplot, file = "pleuro_dotplot.rds")

5 Loading merged pleuro object and ordering the data

Setting the color palette for salamander UMAP, ordering the clusters and generating UMAP and DotPlot.

5.1 Code for general DotPlot

png(filename = "dotplot_fig1.png", width = 10, height = 6, units = "in", res = 400)
DotPlot(pleuro.dotplot, features = c("SNAP25", "SYT1", "RBFOX3", "FOXG1", "SLC17A7", "GAD1","SLC17A6",  "SOX4", "SOX9", "SOX2", "GFAP", "PDGFRA", "NINJ2", "COL1A2","PRSS56", "C1QB", "LCP1"), idents = c("Excitatory telencephalic", "Inhibitory telencephalic", "Non-telencephalic", "Immature neurons", "Ependymoglia", "Oligodendrocyte precursor cells", "Vascular leptomeningeal cells", "Microglia", "Olfactory ensheating cells", "Oligodendrocytes", "Perivascular macrophages")) + RotatedAxis()
dev.off()

null device 
          1 

pdf(file = "dotplot_fig1.pdf", width = 10, height = 6)
DotPlot(pleuro.dotplot, features = c("SNAP25", "SYT1", "RBFOX3", "FOXG1", "SLC17A7", "GAD1","SLC17A6",  "SOX4", "SOX9", "SOX2", "GFAP", "PDGFRA", "NINJ2", "COL1A2","PRSS56", "C1QB", "LCP1"), idents = c("Excitatory telencephalic", "Inhibitory telencephalic", "Non-telencephalic", "Immature neurons", "Ependymoglia", "Oligodendrocyte precursor cells", "Vascular leptomeningeal cells", "Microglia", "Olfactory ensheating cells", "Oligodendrocytes", "Perivascular macrophages")) + RotatedAxis()
dev.off()

null device 
          1 

6 Subsetting Neurons

We subsetted neurons on the basis of neuronal gene expression (SYT1, SNAP25, SLC17A7, SLC17A6, GAD1, GAD2) and re-ran the pipeline for clustering. We RunPCA with 400 pcs, analyzed the ElbowPlot and determined 180 as a good number to calculate UMAP, Neighbors. In terms of clustering, after long deliberation, we settled to res = 6 as it captured most of the cell diversity. Some adjustments to the neuronal object will be made depending on the details of the clustering.

neurons <- subset(pleuro, idents = c('0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '12', '13', '14', '15', '16', '17', '18', '19', '22', '24', '25', '26', '27', '28', '29', '30', '31', '32', '33', '35', '37', '38', '39', '40', '41', '42', '43', '45', '46', '47', '48', '49', '50', '51', '52', '54', '55', '56', '57', '58', '59', '60', '62', '64', '65', '66', '67', '68', '69', '70'))

#neurons <- RunPCA(neurons, npcs = 400)

ElbowPlot(neurons, ndims = 350)

nPCs = 180
  
neurons <- RunUMAP(neurons, dims = 1:nPCs, verbose = FALSE)

neurons <- FindNeighbors(neurons, dims = 1:nPCs, verbose = FALSE)

neurons <- FindClusters(neurons, resolution = 6, verbose = FALSE)

DimPlot(neurons, label = TRUE, reduction = 'umap', label.size = 6, pt.size = 0.7) + NoLegend()

VlnPlot(neurons, features = c("nFeature_RNA", "nCount_RNA"), split.by = "seurat_clusters")
#saveRDS(neurons, file = "neurons_180pcs.rds")

The rest of the neuronal file will be processed in a new R Markdown file.

7 For Glial analysis

We will take only the glial identities in the dataset for a superficial analysis. We are taking all the non-neuronal cells in the dataset and building a heatmap based on known glial gene markers.

oec.markers <- FindMarkers(pleuro.dotplot, ident.1 = "Olfactory ensheating cells")

  |                                                  | 0 % ~calculating  
  |+                                                 | 1 % ~21m 06s      
  |+                                                 | 2 % ~20m 45s      
  |++                                                | 3 % ~20m 31s      
  |++                                                | 4 % ~20m 08s      
  |+++                                               | 5 % ~19m 49s      
  |+++                                               | 6 % ~19m 52s      
  |++++                                              | 7 % ~19m 31s      
  |++++                                              | 8 % ~19m 08s      
  |+++++                                             | 9 % ~18m 47s      
  |+++++                                             | 10% ~18m 27s      
  |++++++                                            | 11% ~18m 09s      
  |++++++                                            | 12% ~17m 51s      
  |+++++++                                           | 13% ~17m 35s      
  |+++++++                                           | 14% ~17m 20s      
  |++++++++                                          | 15% ~17m 04s      
  |++++++++                                          | 16% ~16m 49s      
  |+++++++++                                         | 17% ~16m 35s      
  |+++++++++                                         | 18% ~16m 21s      
  |++++++++++                                        | 19% ~16m 09s      
  |++++++++++                                        | 20% ~15m 57s      
  |+++++++++++                                       | 21% ~15m 44s      
  |+++++++++++                                       | 22% ~15m 31s      
  |++++++++++++                                      | 23% ~15m 21s      
  |++++++++++++                                      | 24% ~15m 10s      
  |+++++++++++++                                     | 25% ~14m 58s      
  |+++++++++++++                                     | 26% ~14m 46s      
  |++++++++++++++                                    | 27% ~14m 34s      
  |++++++++++++++                                    | 28% ~14m 22s      
  |+++++++++++++++                                   | 29% ~14m 09s      
  |+++++++++++++++                                   | 30% ~13m 57s      
  |++++++++++++++++                                  | 31% ~13m 44s      
  |++++++++++++++++                                  | 32% ~13m 31s      
  |+++++++++++++++++                                 | 33% ~13m 19s      
  |+++++++++++++++++                                 | 34% ~13m 08s      
  |++++++++++++++++++                                | 35% ~12m 56s      
  |++++++++++++++++++                                | 36% ~12m 44s      
  |+++++++++++++++++++                               | 37% ~12m 33s      
  |+++++++++++++++++++                               | 38% ~12m 20s      
  |++++++++++++++++++++                              | 39% ~12m 07s      
  |++++++++++++++++++++                              | 40% ~11m 55s      
  |+++++++++++++++++++++                             | 41% ~11m 43s      
  |+++++++++++++++++++++                             | 42% ~11m 30s      
  |++++++++++++++++++++++                            | 43% ~11m 18s      
  |++++++++++++++++++++++                            | 44% ~11m 05s      
  |+++++++++++++++++++++++                           | 45% ~10m 53s      
  |+++++++++++++++++++++++                           | 46% ~10m 41s      
  |++++++++++++++++++++++++                          | 47% ~10m 29s      
  |++++++++++++++++++++++++                          | 48% ~10m 16s      
  |+++++++++++++++++++++++++                         | 49% ~10m 05s      
  |+++++++++++++++++++++++++                         | 50% ~09m 54s      
  |++++++++++++++++++++++++++                        | 51% ~09m 42s      
  |++++++++++++++++++++++++++                        | 52% ~09m 30s      
  |+++++++++++++++++++++++++++                       | 53% ~09m 18s      
  |+++++++++++++++++++++++++++                       | 54% ~09m 06s      
  |++++++++++++++++++++++++++++                      | 55% ~08m 54s      
  |++++++++++++++++++++++++++++                      | 56% ~08m 42s      
  |+++++++++++++++++++++++++++++                     | 57% ~08m 30s      
  |+++++++++++++++++++++++++++++                     | 58% ~08m 18s      
  |++++++++++++++++++++++++++++++                    | 59% ~08m 06s      
  |++++++++++++++++++++++++++++++                    | 60% ~07m 54s      
  |+++++++++++++++++++++++++++++++                   | 61% ~07m 42s      
  |+++++++++++++++++++++++++++++++                   | 62% ~07m 30s      
  |++++++++++++++++++++++++++++++++                  | 63% ~07m 18s      
  |++++++++++++++++++++++++++++++++                  | 64% ~07m 06s      
  |+++++++++++++++++++++++++++++++++                 | 65% ~06m 54s      
  |+++++++++++++++++++++++++++++++++                 | 66% ~06m 42s      
  |++++++++++++++++++++++++++++++++++                | 67% ~06m 31s      
  |++++++++++++++++++++++++++++++++++                | 68% ~06m 19s      
  |+++++++++++++++++++++++++++++++++++               | 69% ~06m 07s      
  |+++++++++++++++++++++++++++++++++++               | 70% ~05m 55s      
  |++++++++++++++++++++++++++++++++++++              | 71% ~05m 43s      
  |++++++++++++++++++++++++++++++++++++              | 72% ~05m 31s      
  |+++++++++++++++++++++++++++++++++++++             | 73% ~05m 19s      
  |+++++++++++++++++++++++++++++++++++++             | 74% ~05m 07s      
  |++++++++++++++++++++++++++++++++++++++            | 75% ~04m 55s      
  |++++++++++++++++++++++++++++++++++++++            | 76% ~04m 43s      
  |+++++++++++++++++++++++++++++++++++++++           | 77% ~04m 31s      
  |+++++++++++++++++++++++++++++++++++++++           | 78% ~04m 20s      
  |++++++++++++++++++++++++++++++++++++++++          | 79% ~04m 08s      
  |++++++++++++++++++++++++++++++++++++++++          | 80% ~03m 56s      
  |+++++++++++++++++++++++++++++++++++++++++         | 81% ~03m 44s      
  |+++++++++++++++++++++++++++++++++++++++++         | 82% ~03m 32s      
  |++++++++++++++++++++++++++++++++++++++++++        | 83% ~03m 20s      
  |++++++++++++++++++++++++++++++++++++++++++        | 84% ~03m 08s      
  |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~02m 57s      
  |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~02m 45s      
  |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~02m 33s      
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~02m 21s      
  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~02m 10s      
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~01m 58s      
  |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~01m 46s      
  |++++++++++++++++++++++++++++++++++++++++++++++    | 92% ~01m 34s      
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~01m 23s      
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~01m 11s      
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~59s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~47s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~35s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~24s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~12s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=19m 37s

8 BarPlot

png(filename = "pleuro_stackedbarplot_animal_Wlegend.png", width = 3, height = 5, units = "in", res = 400)
x
dev.off()

null device 
          1 

pdf(file = "pleuro_stackedbarplot_animal_Wlegend.pdf", width = 3, height = 5)
x
dev.off()

null device 
          1 

pleuro.dotplot <- readRDS(file = "pleuro_dotplot.rds")
data <- table(pleuro.dotplot@meta.data$animal, Idents(pleuro.dotplot))
data_percentage <- apply(data, 2, function(x){x*100/sum(x,na.rm=T)})

library(RColorBrewer)
coul <- brewer.pal(5, "Pastel2") 
ynames <- colnames(data_percentage)

par(oma=c(3,1,1,2) + 0.5)
barplot(data_percentage, 
        col = coul, 
        las =2, 
        font.axis = 2, 
        cex.lab = 0.1, 
        legend.text = rownames(data_percentage),
        args.legend = list(x = "right", inset=c(-0.01,0), xpd = TRUE))

png(filename = "pleuro_barplot_animal_Wlegend.png", width = 30, height = 7, units = "in", res = 400)
par(oma=c(1,1,1,2) + 0.5)
barplot(data_percentage, 
        col = coul, 
        las = 2, 
        font.axis = 2, 
        cex.lab = 0.1, 
        legend.text = rownames(data_percentage),
        args.legend = list(x = "right", inset=c(-0.01,0), xpd = TRUE))
dev.off()

null device 
          1 

pdf(file = "pleuro_barplot_animal_Wlegend.pdf", width = 30, height = 7)
par(oma=c(1,1,1,2) + 0.5)
barplot(data_percentage, 
        col = coul, 
        las =2, 
        font.axis = 2, 
        cex.lab = 0.1, 
        legend.text = rownames(data_percentage),
        args.legend = list(x = "right", inset=c(-0.01,0), xpd = TRUE))
dev.off()

null device 
          1